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Abstract 



We present a new form of explicitly correlated wave function whose parameters are mainly 
linear, to circumvent the problem of the optimization of a large number of non-linear parameters 
usually encountered with basis sets of expUcitly correlated wave functions. With this trial wave 
function we succeeded in minimizing the energy instead of the variance of the local energy, as is 
more common in quantum Monte Carlo methods. We apphed this wave function to the 
calculation of the energies of Be (ls^2p^) and Be" ''S" (ls^2p^) by variational and diffusion 
Monte Carlo methods. The results compare favorably with those obtained by different types of 
explicitly correlated trial wave functions already described in the literature. The energies obtained 
are improved with respect to the best variational ones found in literature, and within one standard 
deviation from the estimated non-relativistic limits. 
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Introduction 



The description of the electron correlation plays a central role in highly accurate quantum 
chemistry calculations. Mean-field methods give a qualitative description for many atomic and 
molecular systems, but in order to get quantitative results the instantaneous correlation between 
electrons must be taken into account. The most common way to include correlation is, starting 
from the Hartree-Fock picture, to approximate the exact wave function using MC-SCF or CI 
expansions. Unfortunately methods based on the orbital approximation converge very slowly to 
the non-relativistic hmit. The reason is that these wave functions include the interelectronic 
distances only in an implicit form. Furthermore this implicit dependence is quadratic instead of 
linear, so the cusp conditions [1] of the exact wave functions are reproduced only for infinite 
expansions. 

A very efficient and effective approach to accurately describe the local behavior of the 
wave function when two electrons colhde is the exphcit inclusion of the interelectronic distances 
into an approximate wave function. Hylleraas [2], Pekeris [3], James and CooUdge [4], and Kolos 
and Wolniewicz [5-7] showed how to obtain very accurate results for two electron systems by 
including the interelectronic distance into the wave function. An alternative possibihty is the 
construction of many-particle permutational symmetry adapted functions in hyperspherical 
coordinates [8, 9]. Unfortunately it is not easy to generalize these methods to many-electron 
systems since the resulting integrals are extremely difficult to evaluate analytically. Beyond four 
electron systems, with at most two nuclei, the analytical approach becomes almost unfeasible [10, 
11]. 

Instead of computing the integrals analytically, one could resort to a numerical method. 
The variational Monte Carlo (VMC) method [12, 13] is a very powerful numerical technique that 
estimates the energy, and all the desired properties, of a given trial wave function without any 
need to analytically compute the matrix elements. For this reason it poses no restriction on the 
functional form of the trial wave function, requiring only the evaluation of the wave function 
value, its gradient and its Laplacian, and these are easily computed. Using the VMC algorithm, 
essentially a stochastic numerical integration scheme, the expectation value of the energy for any 
form of the trial wave function can be estimated by averaging the local energy 

H*I'.p(R)/*I'.p(R) over an ensemble of configurations distributed as , sampled during a 
random walk in the configuration space using MetropoUs [14] or Langevin algorithms [15]. The 
fluctuations of the local energy depend on the quaUty of the function and they are zero if the 
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exact wave function is used (zero variance principle). VMC can also be used to optimize the trial 
wave function and we refer the reader to the hterature for the technical details. 

A popular and effective approach to building compact exphcitly correlated wave functions 
is to multiply a determinantal wave function by a correlation factor, the most commonly used 
being a Jastrow factor [16]. The inclusion of the Jastrow factor does not allow the analytical 
evaluation of the integrals, so the use of VMC is mandatory. However, departing from the usual 
determinantal wave function form can be very fruitful [11], allowing an accurate and, at the same 
time, compact description of atomic and molecular systems. Very few terms are needed to reach a 
good accuracy, in comparison to more common wave function forms. 

The recovery of the remaining correlation energy can be done using the diffusion Monte 
Carlo (DMC) method. Since this method is already well described in the hterature, we refer the 
reader to the available reviews [12, 13]. We only recall here that in this method the exact, but 
unknown, wave function is interpreted as a probabihty density. In the fixed-node (FN) 
approximation [17] the nodal surfaces of the trial wave function are used to partition the 
space and within each region the wave function can be safely interpreted as a probability density. 
It can be shown that the FN-DMC energies are an upper bound to the exact ground state energy. 

This paper is part of an ongoing project in our laboratory to develop accurate and compact 
wave functions for few-electron systems. 

In our previous works [11, 18, 19] we used hnear expansions of exphcitly correlated wave 
functions for calculations on the ground state of few-electron systems. In all cases good VMC 
energies were obtained, both in infinite nuclear mass approximation calculations and non- 
adiabatic calculations. In particular we used a linear expansion of exphcitly correlated 
exponential functions to develop accurate wave functions for two test systems: the beryllium 
atom and the lithium hydride molecule in their ground state. 

Here we present a new form of explicitly correlated wave function and we use VMC to 
extend the apphcation of correlated trial wave function to excited states and five electron systems. 
Furthermore we use DMC to approximate the exact energies and compare them with the 
estimated non-relativistic limits. 

We choose the Be (ls^2p^) and Be" "^8° (ls^2p^) states which are involved in Berylhum 
electron affinity determination. We compare VMC and DMC energies and variances of the 
energy and examine the nodal properties of the trial wave functions comparing FN-DMC results 
with the best variational calculations and the non-relativistic limits estimated by Chung and 
coworkers [20, 21]. 
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Explicitly Correlated functional form 



For an N electron atomic system we write an explicitly correlated trial wave function [11] 



as 



(1) 



In this equation A is the antisymmetrizer operator, (j) is a function of all the electron- 
nucleus distances and g is a function of all the electron-electron distances called correlation 
factor. Both functions include variational parameters. ®s eigenfunction of the spin 

Operators S and of the correct spin multiphcity. The functions (j) and g, being dependent only 

on interparticle distances, are rotationally invariant. This means that their product can describe 
only S states, with zero angular momentum. To describe higher angular momentum states, it is 
necessary to include a function f(r) with the correct rotational symmetry. f(r) is a function of the 
Cartesian electronic coordinates (x,y,z), but might include also the electron-nucleus distances 
[11]. This function might be generahzed including products of the interparticle distances, that 
is ^is the two-body term of a many-body expansion of the wave function. It is possible to further 
generaUze the wave function by taking linear combinations of such terms. 

To assure a high quaUty wave function it is particularly important that the function ^ 
satisfy the cusp conditions [1], representing the behavior of the exact wave function at the 
coalescence of two particles. It is also important to take into account the asymptotic conditions 
[22], which represent the behavior when one of the particles go to infinity. 

The first type of functional form we examined is generated assuming a Pade' factor 

exp[(ar-l-&r^)/(l-l-cr)] for the electron-nucleus part (j) and a Jastrow factor 
exp[a' r/ (1 -I- c' r)] for the interelectronic part g. 



/(r)exp 



ar- + br^ 



exp 



0«. 

S,M, 



(2) 



In the following this wave function will be called Pade'-Jastrow. The Pade' factor is a good 
choice for the electron-nucleus part, because it is the best compromise between flexibihty and 
small number of parameters. In fact this function goes as e'^ for r — > and e^'^'^'^^'for r — > oo . 
So with different exponents it can accommodate both the coalescence at the nucleus and the decay 
for large r. It is also important to point out that this factor can accurately describe both Is and 2s 
orbitals as we have shown in our previous work [11]. 

The main problem with linear expansions of exphcitly correlated trial wave functions is the 
huge number of non-hnear parameters to optimize. In our previous work [11] for more 
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sophisticated factors like Fade' or Jastrow we succeeded in optimizing trial wave functions 
including a maximum of two terms. 

To overcome this problem we choose a second type of functional form, similar to the first 

one: 



/(r)G(r)exp 



+ cr. 









exp 









(3) 



We hmit the expansion to a single term and so we have few non-Unear parameters to 
optimize. However to add extra flexibihty to the wave function we introduce a pre-exponential 
factor G(r) written as a sum of powers of interparticle distances weighted by linear parameters: 

P i q ><j r i j>i s i j 

In the following we will call Eq. 3 pre-exponential wave function (prex). 

Even if this kind of wave function allows us to reduce the effort for the optimization of the 
parameters, we recall that the CPU time needed to evaluate explicitly correlated trial wave 
functions is very large and proportional to the number of permutations generated by the 
antisymmetrizer. No matter the form we choose for exphcitly correlated wave functions, they are 
limited to few electron systems. 

In this paper we also compare Pade'-Jastrow and pre-exponential wave functions with a 
more standard form, widely used in QMC calculations [23, 24], that is the product of a 
determinantal function times the Schmidt-Moskowitz [25] correlation factor (SM): 



(5) 



Det and Det are the determinants for a and P electrons. The function U for atoms is 
given by 



(6) 



where Ck are trial parameters and r = ar/(l+br). 

The determinants are generated from ab initio calculations, in general SCF or MCSCF 
calculations, for a given basis set. Then the correlation factor is added, and its variational 
parameters optimized using VMC calculations. 
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Optimization of the trial wave functions 



Our previous work [11] showed that departing from the usual determinantal wave function 
form can be very fruitful, allowing to write very compact and at the same time very accurate 
wave functions. However it is computationally much more demanding and for this reason special 
care must be given to the design of an efficient way of generating and optimizing the trial wave 
function. These steps must be implemented in the most effective, fast and efficient way. The 
standard way to optimize a trial wave function using VMC is to minimize the variance of the 
local energy using a fixed sample of walkers; a method proposed by Frost [26] and Conroy [27] 
and described in detail by Umrigar, Wilson, and Wilkins [28] and by Mushinsky and Nightingale 
[29]. This has been proved to be numerically much more stable than the energy minimization. For 
our trial wave functions we have found very effective the minimization of the variance of the 
energy 



where the parameter Er can be set equal to the exact energy of the system Eq. Botha(H) and 

|J.(Eo) go to zero as ~^ ' where is the exact eigenfunction: their values for a given 

trial wave function can be used to evaluate the quahty of the trial wave function ^x- 

We used |i(Er) as cost function for the optimization of both Pade'-Jastrow and SCF-MS 
wave functions. 

As to Eq. 2, the optimization of the first term of the expansion is usually performed 
starting from a trial wave function with a reasonable electron-electron Jastrow factor, and with 
the electron-nucleus functions coming from some standard Slater orbital basis set, or from small 
basis sets optimized at the SCF level. As we showed in our previous work [11], it is possible to 
build a trial wave function as a linear expansion of n terms by adding an extra term to an 
optimized n-1 term wave function. This procedure worked well for two and three electron systems 
with simple exponential basis sets, but not for more sophisticated Fade' and Jastrow basis sets. 
For these reasons in this work we optimized only one term Pade'-Jastrow functions. 

Let us now consider the case of the pre-exponential trial wave function. The function in 
Eq. 3 can be written as a hnear combination: 




(7) 



or, even better, of the second moment with respect to an arbitrary parameter Er, |i(Er): 




(8) 



(9) 
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where the term is given by 



(10) 



<I>(p) is a single term Pade'-Jastrow function whose parameters p are optimized minimizing 
|a.(ER) in a preliminary step, and then are fixed during the optimization of the linear parameters. 

As to the linear parameters, we succeeded in minimizing the energy instead of the variance 
of the local energy or the second moment |J-(Er), so we could choose the best linear parameters 

according to the observable we are interested in. The standard linear variational methods requires 
the solution of the secular problem, and so the calculation of the matrix elements 



These integrals are evaluated during a VMC simulation. 

We show the main features of this pre-exponential trial wave function using as benchmarks 
the Be and Li ground states. 

In table 1 we compare the energy for Be 'So ground state obtained by the one term Pade'- 
Jastrow function with the value calculated by the function obtained adding a pre-exponential 
including all the electron-nucleus and electron-electron distances and their products, a total of 66 
terms. The two values evidence a large gain of correlation energy for the addition of the pre- 
exponential factor and further improvement is obtained adding 20 more terms, that is the third 
and fourth powers of the electron-nucleus and electron-electron distances, to the linear expansion. 
So the pre-exponential factor adds flexibihty to the wave function in a very efficient way. 

To examine the relative efficiency of linear and non linear parameters in adding flexibility 
to the wave function, for the Li ground state (see table 2) we compare the energy of a 28-term 
pre-exponential function with the result of a 8 -term expansion of explicitly correlated exponential 
functions, a simplified form of the Pade'-Jastrow wave function in which both the factor (|) and the 
correlation factor g are in the form exp(cr). These two trial wave functions give a similar gain of 
correlation energy, but in the pre-exponential case there are 12 non-hnear and 28 linear 
parameters, while the 8 term expansion of correlated exponential functions includes 48 non-hnear 
and 8 linear parameters. In spite of the smaller number of parameters the pre-exponential 
function gives a better result, and its optimization process was much easier and faster. 

We optimized the linear parameters also by minimization of |J-(Er). The calculated energy 
at VMC level is worse than the one obtained by minimization of the energy: this is obviously 
related to the different minima of the energy and |i(Er). 



(11) 




8 



Results and discussion 



For Be (ls^2p^) and Be" ^S" (ls^2p^) we computed SCF and CASSCF trial wave 
functions using GAMESS with the Slater orbital basis sets reported in table 3. Each orbital was 
fitted with 6 Gaussian fiinctions. We optimized the Schmidt-Moskowitz correlation factor of the 
SCF-SM functions and the non-linear parameters of the Pade'-Jastrow and pre-exponential 
functions minimizing the variance of the local energy. 

Beside the VMC energies, we report the variance of the local energy, given by Eq. 7, 
estimated using VMC. DMC energies were obtained by a linear fit of the energy at three time 
steps (x=5, 3, and 1 mhartree"^) and extrapolation to x=0 mhartree"\ 

We compare our results with the best variational energies obtained by Chung and 
coworkers [20, 21], who used linear expansions of Slater orbitals in the L-S'-couphng scheme. 

Be (ls^2p^) 

The calculations for the excited state Be ^P (ls^2p^) were carried out with the following 
three trial wave functions: 

1) a single determinantal function times a nine term Schmidt-Moskowitz correlation factor; 

2) a one term Pade'-Jastrow function; 

3) a pre-exponential function with 33 terms. 

The spin eigenfunction used for the Pade'-Jastrow and pre-exponential function is 
0f 1 = {ocfiococ — ficxCXCX) . The pre-factor f(r) that defines the state symmetry is 

f(r) = x^y^-x^yj (12) 
The results are reported in table 4. 

As to the determinantal wave function, we used only a single determinant, as the energy 
lowers by only 0.003 hartree on going from the SCF to a CASSCF function for two electrons in 
an active space of two p shells (15 configurations, the highest weight of the first double excitation 
being equal to 0.04). 

The exphcit inclusion of the interelectronic distances in the wave function by the SM factor 
results in a large improvement of the quahty of the wave function, as shown from the lowering of 
the energy (0.043 hartree) and the variance of the energy (almost four times smaller) on going 
fi-om SCF to SCF-SMA^C. 
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A further improvement of the trial wave function at VMC level is found using more 
sophisticated functional forms [11], like the Pade'-Jastrow and the pre-exponential ones. For the 
pre-exponential function we used a 33 term expansion, including all ri, ti, and the products riry. 

In particular in the case of the pre-exponential function we were able to optimize the hnear 
parameters of G(r) minimizing the energy, not the variance of the local energy, obtaining an 
energy 0.9 mhartree higher than the best variational one. We also notice that this wave function is 
very compact with 45 (12 non-linear and 33 linear) variational parameters on the whole. 

At DMC level already the SCF-SM wave function gives a lower energy than the best 
variational value, and 0.2 mhartree higher than the estimated non-relativistic hmit (NRL). It 
means that the nodal surfaces of this function are fairly good, at variance with the SCF-SM trial 
wave function for the Be ground state, whose energy is 11 mhartree higher than the NRL [30]. 
This large nodal error is due to the strong contribution of the first double excitation in improving 
the quahty of the nodal surfaces, because of the quasi-degeneracy of the 2s and 2p orbital. 

DMC energies for Pade'-Jastrow and pre-exponential functions have the estimated NRL 
within one standard deviation, that is the nodal surfaces of these wave functions are correct and 
better than the SCF ones. 

Be V (ls^2p^) 

The calculations for the excited state Be" "^S" were carried out with these four trial wave 
functions: 

1) a single determinantal function times a nine term Moskowitz-Schmidt correlation factor; 

2) a multideterminantal function times a nine term Moskowitz-Schmidt correlation factor; 

3) a one term Pade'-Jastrow function; 

4) a pre-exponential function G(r) with 61 terms. 

The CASSCF wave function for three electrons in an active space of two p shells includes 
20 configurations. The first two highest weights, relative to the first double and single excited 
configurations, are equal to 0.125 and 0.03 and indicate a more marked multiconfigurational 
character of the wave function. 

As we have seen for Be state, the gain in energy and the lowering of the variance of the 
energy between SCF and SCF-SM are very large, while between SCF-SM and CASSCF-SM 
they are an order of magnitude less. 

The spin eigenfunction used for the Pade'-Jastrow and pre-exponential function is 
©3/23/2 = (CcPaaa - Paaaa) . The pre-factor f(r) that defines the state symmetry is 

f (r) = X3y4Z5 + X5y3Z4 + X4y5Z3 - X3y5Z4 - X5y4Z3 - X4y3Z5 (13) 
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The results are reported is table 5. 

For Be" using the Fade' and pre-exponential functions we obtained better energies and 
variances of the energy than the SCF-SM and CASSCF-SM ones, as already seen in the Be 
case. In particular, with the pre-exponential function with 61 linear parameters G(r) (all ri, ti 
and the products rir^), the VMC energy is 0.9 mhartree higher than the best variational energy. 

From DMC simulations we see clearly that the nodal surfaces don't change on going from 
the SCF to the CASSCF trial wave function and in both cases we have around 0.4 mhartree of 
nodal error. 

Pade'-Jastrow and pre-exponential functions have better nodal surfaces and their DMC 
energies have the estimated NRL within one standard deviation. 

Conclusions 

We have used exphcitly correlated functional forms to improve the quahty of the trial wave 
functions usually adopted to calculate the energy of a system. For the two excited states of Be 
and Be" we obtained better non-relativistic energies with very compact trial wave functions 
compared to the best variational results. 

Using a suitable pre-exponential factor we were able to improve the flexibihty of the trial 
wave function without including too many non-hnear parameters: this kind of trial wave function 
allowed us to minimize directly the energy instead of the variance of the local energy. As to the 
computational time, the optimization of the Be five electron trial wave function and the VMC 
calculation required around a week on a modern PC .It is not possible to compare this CPU time 
with calculations by correlated Gaussians as at present they are limited to four electron systems. 

Our DMC energies are in good agreement with the estimated NRL obtained by Chung and 
coworkers [20, 21]. 

From our best values for Be ^P (-14.39547(5) hartree) and Be" (-14.40620(6) hartree) 
we compute an electron affinity of 0.01073(8) hartree = 292(2) meV, within two standard 
deviations from the experimental value 295.49(25) meV [31] and the theoretical value 295.0(7) 
calculated by Hsu and Chung [20]. A significant comparison would require the reduction of the 
calculated standard deviation by one order of magnitude. 
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Table 1 . VMC results for Be ground state. 





VMC Energy (hartree) 


% Correlation energy 


one term Pade'-Jastrow 


-14.6528(2) 


84.57 


66 term prex 


-14.6633(3) 


95.70 


86 term prex 


-14.6651(2) 


97.60 


HF limit 


-14.57302 




NR limit 


-14.66735 
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Table 2. VMC results for Li ground state. 



VMC Energy (hartree) % Correlation energy 



Stermexp. -7.4775(2) 98.29 

28termprex -7.47770(8) 99.20 

HF limit -7.43274 

NR limit -7.47806 
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Table 3. Basis sets for SCF and CASSCF calculations 



System 


r$e r 


T3„- 4c;o 

rie O 


1 n 

Is 


D. 1 


J. / 


1 ^ 

Is 


4.2 


4.2 


2s 


4.3 


4.3 


2s 


2.4 


2.4 


2p 


1.65 


1.65 


2p 


0.76 


0.76 


2p 




0.376 
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Table 4. Be (ls^2p^) energies and VMC variances of the energy 



Method 


Energy (hartree) 


<'vMc(H) 


SCF 


-14.3340 


1.68(2) 


SCF-SMAHVIC 


-14.3769(2) 


0.48(1) 


Pade'AHVlC 


-14.3930(1) 


0.27(1) 


33 term prex /VMC 


-14.3942(1) 


0.22(1) 


Best variational 


-14.3951086 




SCF-SM/DMC 


-14.39521(5) 




Pade'/DMC 


-14.39541(7) 




33 term prex /DMC 


-14.39547(5) 




Estimated LNR 


-14.3954404 
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Table 5: Be ^S" (ls^2p2) energies and VMC variances of the energy 



Method 


Energy (hartree) 


<'vMc(H) 


SCF 


-14.326976 


1.68(1) 


CASSCF 


-14.334010 


1.68(1) 


SCF-SMAHVIC 


-14.3769(2) 


0.48(1) 


CASSCF-AHVIC 


-14.3836(1) 


0.48(1) 


Pade'AHVlC 


-14.4031(2) 


0.29(1) 


61 term prexAHVIC 


-14.4051(2) 


0.21(1) 


SCF-SM/DMC 


-14.40594(8) 




CASSCF-SM/DMC 


-14.40597(7) 




Best variational 


-14.4060320 




Pade'/DMC 


-14.40620(4) 




61 term prex /DMC 


-14.40620(6) 




Estimated LNR 


-14.406282(26) 
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